function eq = solve_eq_ss(param,glob,options)

% Finds initial steady state
% With consumption normalized to 1

 D1 = 0;
 D2 = 50;
 Dmax = D2;
 
% Impose that C = 1. 
param.C = 1;   

Ds    = nan(options.itermax_DC,1);

if nargin == 5
    param.D = eq0.D_out;
    options.cresult = eq0.c;
else
    param.D = .459;
end
    
    
for t = 1:options.itermax_DC
    
    param.D = (D1 + D2)/2;
    Ds(t) = param.D;        

    param.w = param.omega*param.C; % HH FOC

    tic
    eq = solve_cL(param,glob,options);     % solve incumbent & entrant pbms  
    eq = setup_transition_matrices(eq, param,glob,options); % set up transition matrices
    [diff,eq] = find_stationary_dist(eq, param,glob,options);
    toc
        
   if abs(diff) < options.tolD
        break
   end
   
   disp(diff)
   
   D1s(t) = D1;
   D2s(t) = D2;
   Ds(t) = param.D;
   
   subplot(3,1,1)
   hold off
    plot(1:t, D1s, 'rx', 'MarkerSize', 5); drawnow;
    hold on
    plot(1:t, D2s, 'rx', 'MarkerSize', 5); ylim([0 Dmax]); drawnow;
    plot(Ds); title('Guessed value of demand index'); drawnow;

    subplot(3,1,2)
    diffs(t) = diff;
    plot(1:t, diffs); title('Expeted value of entry'); drawnow;
    
   if diff < 0
       D2 = param.D; 
    else
       D1 = param.D;
   end
   

    
end

    subplot(3,1,3)
    plot(glob.bgridf, eq.Lp)
    title('Distribution of labor'); drawnow;
    options.cresult = eq.c;
 
   
%    param.D = options.damp_DC*param.D + (1-options.damp_DC)*eq.D_out;
    eq.D = param.D;
    
    eq.flag = 0;
    if t == options.itermax_DC
        eq.flag = 1;
    end
    
end